perm filename SSS.F4[MUS,LCS]1 blob
sn#007339 filedate 1974-01-08 generic text, type T, neo UTF8
SUBROUTINE SSS(VV,N1,A1)
DIMENSION V(50,4),A1(512),C(30,4),YP(30),J(30),NX(3),KA(14),K(9)
DIMENSION VV(50,4)
EQUIVALENCE(K1,K(1)),(K2,K(2)),(K3,K(3)),(K4,K(4)),(K5,K(5)),
1 (K6,K(6)),(K7,K(7)),(K8,K(8)),(K9,K(9))
DATA KA/1,2,2,1,1,2,1,1,0,2,1,-1,0,1/,DX/.00001/
IF(VV(1,2).EQ.0) VV(1,2)=1
DO 5 I=1,30
DO 5 L=1,2
5 V(I,L)=VV(I,L)
NX(1)=N1
698 NX(2)=NX(1)-1
DO 10 I=1,NX(1)
10 V(I,2)=(V(I,2)-1)/99.
DO 20 I=2,NX(2)
YP(I)=(V(I+1,1)-V(I-1,1))/(V(I+1,2)-V(I-1,2))
20 IF((V(I+1,1)-V(I,1))*(V(I,1)-V(I-1,1)).LE.0) YP(I)=0
DO 22 I=1,9
22 K(I)=KA(I)
KOUNT=0
21 KOUNT=KOUNT+1
V1=V(K2,1)-V(K1,1)
V2=V(K2,2)-V(K1,2)
802 IF((YP(K2)-V1/V2)*(V(K3,1)-V(K4,1)).GT.0) GO TO 30
24 Z=V(K2,K5)+(V(K1,K6)-V(K2,K6))*YP(K2)**K7
IF(YP(K2)**2.LT.DX.AND.V1**2.LT.DX) GO TO 36
IF(YP(K2)**2.LT.DX) GO TO 38
D1=V(K2,K5)-Z
806 D2=Z-V(K1,K5)
ZZ=(V(K1,K6)*D2+V(K2,K6)*D1)/(D1+D2)
808 YP(K1)=(ZZ*K9+V(K2,1)*K8-V(K1,1))/
1 (ZZ*K8+V(K2,2)*K9-V(K1,2))
GO TO 40
30 DO 32 I=5,9
32 K(I)=KA(I+5)
GO TO 24
36 YP(K1)=0
GO TO 40
38 YP(K1)=-100
IF(KOUNT.EQ.2) GO TO 39
IF(V(K2,1).GT.V(K1,1)) YP(K1)=100
GO TO 40
39 IF(V(K2,1).LT.V(K1,1)) YP(K1)=100
40 IF(KOUNT.EQ.2) GO TO 50
DO 42 I=1,2
K(I)=NX(I)
42 K(I+2)=K(I)
DO 44 I=5,9
44 K(I)=KA(I)
GO TO 21
50 NX(3)=NX(2)-1
N=1
52 N=N+1
IF(N.GT.NX(3)) GO TO 92
V1=V(N+1,1)-V(N,1)
V2=V(N+1,2)-V(N,2)
Y1=YP(N)-YP(N+1)
IF(Y1**2.LT.DX.AND.V1**2.GT.DX) GO TO 720
710 X=(V1-YP(N+1)*V(N+1,2)+YP(N)*V(N,2))/Y1
715 IF(X.GE.V(N,2).AND.X.LE.V(N+1,2)) GO TO 52
IF(Y1**2.LT.DX.AND.V1**2.LT.DX) GO TO 52
720 DO 120 I=NX(1),N+1,-1
V(I+1,2)=V(I,2)
V(I+1,1)=V(I,1)
120 YP(I+1)=YP(I)
YP(N+1)=.5*V1/V2
IF(V1*(YP(N)-V1/V2).LE.0) YP(N+1)=4*YP(N+1)
V(N+1,2)=.5*(V(N+2,2)+V(N,2))
V(N+1,1)=.5*(V(N+2,1)+V(N,1))
N=N+1
DO 88 L=1,3
88 NX(L)=NX(L)+1
GO TO 52
92 DO 140 I=1,NX(2)
W0=YP(I)
W1=YP(I+1)
W2=V(I+1,2)-V(I,2)
W3=V(I+1,1)-V(I,1)
C(I,1)=(W2*(W0+W1)-2*W3)/(W0-W1)
C(I,2)=W2-C(I,1)
C(I,4)=W0*C(I,2)
140 C(I,3)=-C(I,4)+W3
730 DO 150 I=1,NX(1)
150 J(I)=511*V(I,2)+1
740 DO 160 I=1,NX(2)
L1=J(I)+1
IF(I.EQ.1) L1=1
L2=J(I+1)
750 DO 160 L=L1,L2
X=(FLOAT(L)-1.)/511.
IF(C(I,1)**2.LT.DX) GO TO 155
ZX=.5*SQRT(C(I,2)**2-4*C(I,1)*(V(I,2)-X))/C(I,1)
T1=-.5*C(I,2)/C(I,1)+ZX
T2=T1-2*ZX
IF(T2.GT.-DX.AND.T2.LT.(1+DX)) T1=T2
155 IF(C(I,1)**2.LT.DX) T1=-(V(I,2)-X)/C(I,2)
160 A1(L)=C(I,3)*T1**2+C(I,4)*T1+V(I,1)
770 END